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The leading-order nucleon-nucleon (NN) potential derived from chiral perturbation theory con- 
sists of one-pion exchange plus short-distance contact interactions. We show that in the 1 So and 
3 Si- 3 Di channels renormalization of the Lippmann-Schwinger equation for this potential can be 
achieved by performing one subtraction. This subtraction requires as its only input knowledge of 
the NN scattering lengths. This procedure leads to a set of integral equations for the partial-wave 
NN t-matrix which give cutoff-independent results for the corresponding NN phase shifts. This 
reformulation of the NN scattering equation offers practical advantages, because only observable 
quantities appear in the integral equation. The scattering equation may then be analytically con- 
tinued to negative energies, where information on bound-state energies and wave functions can be 
r*j ' extracted. 
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O 1 I. INTRODUCTION 

Chiral perturbation theory (xPT) G], HI 9 provides a systematic technique by which nuclear forces can be derived. 
As first pointed out by Weinberg [4|, Q , the nucleon-nucleon (NN) potential can be expanded in powers of the xPT 
' expansion parameter P = P, J^"" . The leading-order (LO=P°) potential is then the venerable one-pion exchange (OPE), 

■ supplemented by a short-distance "contact interaction" , which is a constant in momentum space, or equivalently in 
co-ordinate space a three-dimensional delta function. Thus, the LO xPT potential is singular, and when iterated in 

. a Lippmann-Schwinger (LS) equation requires regularization and renormalization. 
t— I ' Weinberg's idea has since been extended to several orders in the chiral expansion with NN potentials of x?T derived 

! to NNLO in Refs. @, 0, and to N 3 LO @, At n3l0 a fit to NN data bclow T iab = 200 McV can be achieved 
which is comparable in quality to that obtained in 'high-quality' NN potential models, provided cutoffs in the range 
f"-- | 500-600 MeV are considered [fj. For this cutoff range the phase shifts obtained from the N 3 LO xPT potential have 

■ little residual cutoff dependence at Ti a b < 200 MeV. 
However, cutoff independence over a wider range is desirable. In field-theoretic terms the issue of whether this can 

be achieved or not speaks to whether a single momentum-independent contact interaction is sufficient to renormalize 
the singular OPE potential [H|,[l3|- Beane, Bedaque, Savage, and van Kolck imposed a co-ordinate space regulator on 
the LO xPT potential, and adjusted the two free parameters that multiply the contact interactions so as to reproduce 
both S-wave NN scattering lengths. They then showed that the resulting phase shifts are independent of the regulator 
scale (up to small higher-order corrections) when it is varied over an order of magnitude [l3|. (Ref. [l3[ also showed 
that renormalization of the m^-dependence of the nuclear force requires the contact interaction to be a function of 
m w , but such m^-dependence will not concern us here.) This cutoff-independence of NN S-wave phase shifts resulting 
from the LO potential derived by Weinberg has since been confirmed by Ruiz Arriola and Pavon Valderrama [3], 
and — using a momentum-space cutoff — by Nogga, Timmermans and van Kolck [l5[. In addition, Ref. [l5[ examined 
the behavior of this LO potential in higher partial waves and argued that channels with I > 1 in which the OPE 
potential is attractive also need additional contact interactions if they are to be properly renormalizcd, a finding 
which was confirmed in Ref. (l6| . This analysis agrees with an examination of the renormalization-group scalin g of 
short-distance operators in the presence of one-pion exchange 17 1. The usefulness of the results of Refs. [IH [r^,I17T| 



has, however, been disputed in Ref. [18(. For a review of the situation see Ref. [1| 
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However, the NN partial waves on which this ongoing discussion centers are not the focus of this study. Instead, 
we will re-examine the renormalization of the leading-order %PT potential in the 1 So and 3 Si- 3 D! channels. We will 
develop a strategy that replaces the divergent integrals and unphysical contact interactions that appear in the LS 
equation for this potential with finite integrals and physical quantities. In particular, we will show that knowledge of 
the zero-energy on-shell NN scattering amplitude allows us to turn an ill-behaved LS equation into a set of equations 
that give cutoff-independent results which are equivalent to the ones obtained in Ref. [la ]. 

Such a conversion of divergent loop integrals and potentially infinite Lagrangian couplings into a sum of experimental 
quantities and convergent integrals is a well-known strategy in perturbation theory. In a renormalizablc theory it was 
shown to be possible to pursue this strategy to all orders in perturbation theory by Bogoluibov, Parasuik, Hepp, and 
Zimmerman 2fjj . Hereafter we refer to this approach as "subtractive" renormalization. We show that such subtractive 
renormalization can be applied to the LS equation for the LO xPT potential. This is a non-trivial result, as while 
contact interactions have simple integral representations in perturbation theory the contact terms appearing in the LO 
xPT potential do not. In consequence the development of subtractive renormalization in a non-perturbative context 
is, to our knowledge, only a recent advance, having been originally developed for the "pionless EFT" in Refs. U 221. 

An earlier attempt to develop subtractive renormalization for the NN interaction of ;\PT was made in Ref. 231 ] . 

and extended in Ref. [24[. There a large negative energy, E sub = —jj, is used as subtraction point while solving 
the LS equation at a fixed cm. energy in the NN system. The assumption that the Born approximation holds at 
E su by provided the parameter fi is large enough, then determines the behavior of the NN amplitude at that energy. 
However, this assumption is not generally valid for a singular potential. Frederico et al. therefore effectively enforce 
a particular off-shell dependence of the t-matrix at E su i,. We will show that such an assumption is not necessary. 
Instead, subtractions first suggested by Hammer and Mehen [2l| can be used at E = to obtain the half-shell NN 
t-matrix t(p, 0; E = 0) from the original LS equation and the NN scattering length. From this, NN phase shifts at any 
energy can be be obtained assuming o nly Hermiticity of the underlying potential. The latter technique was developed 
and employed by Afnan and Phillips [22] in the context of the low-energy three-body problem, but was not yet used 
in a two-body context. 

The paper is structured as follows. In Section [TT] we review the LO xPT NN potential and solve the LS equation 
for S-waves employing a momentum cutoff. We demonstrate the well-known result that S-wavc phase shifts are cutoff 
dependent unless the contact interaction is adjusted so as to reproduce some low-energy observable. We reproduce the 
results of Ref. [l5| for cutoffs in the range 0.5 to 4 GeV. In Sec. IHII we perform subtractions on the LS equation such 
that reference to the contact interaction disappears and demonstrate that this procedure results in cutoff-independent 
predictions for NN phase shifts. Our use of subtractive renormalization allows us to consider cutoffs that are much 
larger than those used in Ref. [l5[, and in Sec. Mil we consider A as high as 50 GeV in order to show what is possible 
with our new technique. We also explain the differences between our approach and that of Ref. [23j|. In Sec. IIVI we 
analytically continue the subtracted equation to negative energies and obtain bound-state energies and wave functions, 
thereby showing that our techniques arc equally suitable for obtaining the NN bound state energy and wave function. 
We summarize and outline possible extensions of this work in Sec. [V] 



II. THE LEADING-ORDER NN POTENTIAL IN CHIRAL PERTURBATION THEORY 

As first pointed out by Weinberg in Refs. [3, H[ , in x?T t ne operators of lowest chiral dimension which can contribute 
to the NN potential are: 

(p'\V\p) = ~n ■ ^ ai '% > ~_ P py + %~ P) + '[(Or + 3C S ) + n ■ r 2 (C s - C T )} (2.1) 

where ti, t-i are the isospin and ci, o-i the spin operators. The couplings constants are given by gA = I -25 and / ff = 93 
MeV, the cm. momentum by p. The low-energy constants Cs and Ct represent the effects of degrees of freedom 
with masses of order A x . Their value is not determined by chiral symmetry, and they must be fitted to experimental 
or lattice [25| data. Weinberg then advocated inserting this potential into the non-relativistic Lippmann-Schwinger 
equation, which is sufficient because we arc considering a low-energy expansion of V: 

f d 3 v" 1 
(p'\t(E)\p) = (p'\V\p) + / P(p'\ V \p")— -^(p'^^lp) (2.2) 

(with e a positive infinitesimal and M the nucleon mass). This resums the effects of infra-red enhancement due to 
the presence of heavy (cf. m n ) particles in the intermediate NN state. In a partial-wave expansion Eq. (|2.2[) reads 

tvi(p,p;E) = wi(p ,p) + y.-M j— — , (2.3) 

„ TT J a p^ + ie- p' u 
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where Pq/M = E is the cm. energy, and I denotes the angular-momentum quantum number. For the singlet state 
we have I = V = I" = 0, for the triplet state these values can be or 2. The explicit expressions for the partial wave 
projections of OPE are given in Ref. [23| . For the convenience of the reader we repeat them here. For the 1 So state 
one obtains (with Ct,Cs = 0) 



2 / A 
9a 



^"'igl'-j- v^-;,.^ ]' (2 - 4) 

The angular- momentum projected matrix elements for the 3 Si- 3 Di channels (1,1' = 0,2) are given by 



327r/2 y J_ x p>2 +p 2 — 2p'px + m 2 
2 V2 / f 1 , p' 2 -2p'px+p 2 ax 2 - i) 



v * (p' t p) = m± \ dx P -'PPX+P W -2 

IGttJ 2 p' 2 + p 2 — 2p'px + m 2 

t , 9l ( f 1 . 2ppix-{p 2 +p' 2 )(±x 2 -\) \ 

V 2 2{P,P) = oo (2 / dx 72~, 2 ^1 i 2 ■ ( 2 - 5 

327r/^ \J_ 1 p' 2 + p 2 - 2p'px + mi 



Inspection of the limits for p,p' — > oo for v^(p' ,p) in Eq. (|2.4|) gives [2S 



hm _y m tf,p) = J!L s . (2.6) 



p or p'— too lOTT j. 



For the triplet channels the following limits give a finite value [22 
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lim t4(p',p) 



167T/2' 

vl Q {p',p) = 9 -^§, (2.7) 

while linip—too v| (p',p) = hmp.p'-^oo v\ 2 ip \p) = 0. Thus, the integrals of Eq. (|2.3p diverge. The ultra-violet behavior 
of the integral equation can be regulated by introducing a momentum cutoff A as the upper bound of the integration. 
Meanwhile, the infra-red singularity at p = po in Eq. (|2.3[) is separated into a principal-value part and a delta-function 
part, with the principal- value singularity treated by standard subtraction techniques [26| . 

Once these two steps are made the integral equation can be solved by standard Gauss-Legendre quadrature. The 
results of our numerical calculations for the 1 So and 3 Si- 3 Di phase shifts are shown on the left-hand side of Figs.Q] 
. We confirm that if Cs{Ct) takes a fixed value (e.g. zero), the resulting 1 Sq ( 3 Si) phase shift has sizable cutoff 
variation — of order <5i5 {63$^ itself — even for small projectile laboratory energies Ti a b. 

One way to remedy the situation is to renormalize each of the S-wave potentials by taking advantage of the presence 
of the arbitrary constants Cs and Ct- In Ref. [l5[ Eq. (|2.2[) was solved with the relevant constant adjusted such that 
for a given cutoff A the S-wave phase shifts at 10 MeV were reproduced. Here we do something similar, except that 
we use the NN scattering lengths as input, adopting the values a s = —23.7 fm for the singlet case and at = 5.43 fm 
for the triplet [13 . 

Once Cs and Ct are fitted in this way we can predict the phase shift at an arbitrary energy, and the resulting 
"renormalizcd" phase shifts are plotted in in the right-hand panels of Figs. Q]as a function of the laboratory kinetic 
energy. The phase shifts are now approximately cutoff independent, with the dependence on A relegated to effects 
~ 1/A in the effective ranges, r s and r t (see Table [J). This reproduces results of Beane, Bedaque, Savage and van 
Kolck and others 0EH]. 

The renormalization constants Cs and Ct that were used to generate these results are given in Table |TT] together 
with the corresponding cutoff A. While Cs exhibits a monotonic behavior as function of A, the evolution of Ct with 
A is particularly interesting, and is shown in Fig. [2] This confirms the result shown in Fig. 2 of Ref. [l5j|, and shows 
that Ct(A) diverges at specific values of A. The origin of this behavior is the appearance of additional bound states 
in the 3 Si- 3 Di channel as the cutoff is increased. At the specific values of A where one more bound state appears, 
Eq. (|2.2p cannot be solved by fitting Cr(A) to reproduce a t . 

But even when the constants Ct(S) arG finite, fitting them in order to reproduce the scattering length can be very 
delicate. For example, at A = 1 GeV Cs must be tuned to 1 part in 10 4 in order to get a s accurate to 1%, while 
for A = 10 GeV it must be tuned to 1 part in 10 6 to achieve this accuracy. At arbitrarily large cutoffs it therefore 
becomes numerically impossible to fit Cs and Ct- 



4 



III. SUBTRACTED INTEGRAL EQUATIONS 

An alternative to fitting the contact interactions C5 and Ct so as to reproduce the 1 So and 3 Si scattering lengths 
is to perform sub-tractive renormalization on the LS equations for the potential (|2.f [) . In this section we explain the 
subtractive renormalization procedure. Note that in contrast to the approach advocated in Ref. (23[, our subtraction 
procedure does not require assumptions about the behavior of t anywhere. 

Let us first concentrate on the cm. energy E = in the NN system. The definition of the scattering lengths a s M 
in the limit E = p 2 /M — > gives in the singlet (triplet) states 



t s{t) (0 E - 01 - } - ^1 
t 0Q [U,U,E-U)- _ MpQ - M , 



where foo\E) represents the scattering amplitude in the singlet (s) or triplet (t) channel. The half-shell and on-shell 
LS equations for this energy read 

t{p, 0; 0) = v(p, 0) + C + £ dp' p' 2 ( " (P '^ /2 +C ) t(p', 0; 0) (3.2) 
i(0, 0; 0) = v(0, 0) + C+ —M [ A dp' p' 2 f v (°^) + C \ t(p / ) . 0)) (3 .3) 



where we dropped the channel indices since our procedure is valid for the singlet as well as triplet LS equations. 
Subtracting Eq. (j33]) from Eq. (|3~2")) and using (f3~i"j) leads to 



in the singlet channel, while for the triplet channel the equation corresponding to Eq. (|3.4p is 

Up', 0; 0) - J2 hl£ dp"p" 2 ( 4 '" (P>, ^ /2 4 '" (0,p/,) ) Uip", 0; 0) = vj n (p', 0) - v^O, 0) + & m 5 m ± (3.5) 

where 1,1' = 0,2. Note that since «2o(0,p;0) = v^ip, 0; 0) = 0, and 1^2(0, P',0) = «22(p, 0; 0) = 0, the corresponding 
t-matrix elements £02 and £22 must also be zero. 

Equations akin to Eqs. (|3.4p and (|3.5p were first developed for the three-body problem in Ref. [2l|. They have the 
advantage that the constants Cs and Ct cancel out in the differences that form Eqs. (|3.4[) and (|3.5|) . so that these 
equations contain only observable quantities. They can be solved using standard techniques. Our results shows the 
relative difference between the two methods is less than 0.2% for all channels. 

Having determined t(p, 0; 0) we set about obtaining half-shell t-matrix elements i(p,p';0), p' being an arbitrary 
momentum, by carrying out a further subtraction. First we note that: 

t( P ,p';0) = V (p,p>) + C+*M f A dp" P " 2 ( V{P,P "l + C ) tip",p';0) (3.6) 



k Jo V -P" 

rA 



t(0,p';0) = v(0,p') + C+hl J dp"p" 2 ( V{0, ^ C y( P ",p';0). (3.7) 
Subtracting Eq. (f3~7|> from Eq. (|3~6)l leads to 

tip,p';Q)-ti0,p';0)=vip,p')-vi0,p') + ~M £ dp" p" 2 MEi^I ~^ P " ] ^ t(p",p'; 0). (3.8) 

Taking advantage of the property tip', 0; 0) = i(0,p'; 0) T , we can use tip', 0; 0) obtained from Eq. (|3.5|) to solve for 
tip,p'; 0) in Eq. (|3.8p . The calculated results are displayed in Fig.[3]and show good agreement (relative difference^ 2%) 
between this subtractive method and the one described in the previous section. 

As a last step we need to obtain the general half-shell t-matrix tip,p'; E) at arbitrary cm. energies E from the one 
at zero energy that is given by Eq. (|3.8[) . We start again from two LS equations, 

tip,p';0) = vip,p') + C+^M [dp"p" 2 ( " (P '^L +C ) tip', P ";0). (3.9) 
tip,p';E) = vip,p')+C+-M f A dp" P " 2 ( V ^[l + l C . ) t(p',p";E), (3.10) 

n JO \PQ ~ P + * e / 
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which, in order to express the idea more clearly, we simplify to 



f(0) = (v + C) + (v + C)g(Q)t(Q), 
t(E) = (v + C) + (v + C)g{E)t{E). 



The second equation, Eq. (|3.12p , can be rewritten as 

t(E)(l +g(E)t(E))- 1 =v + C, 
whereas the first, Eq. (|3. 1 1[) . can also be expressed as 

f(0) = (v + C)+t(0)g(0)(v + C). 
Substituting Eq. (|3~T3l) into Eq. (|5T3|t gives 

i(0) = [1 + t(0) fl (0)] t(E) [1 + giE^E)}- 1 

which after simplification yields 

t(E) = t(0)+t(0)(g(E)-g(0))t(E). 



(3.11) 
(3.12) 



(3.13) 
(3.14) 
(3.15) 
(3.16) 



We can now obtain t(E) from t(0) by solving the integral equation, Eq. (|3.16|) . Then choosing p = po with E = Pq/M 
and setting p' = p we have the on-shell t- matrix element iz(poi Po! E) at an arbitrary energy E. The on-shell element 
is related to the phase shift Si through the well-known relation 



ti{p ,p ;E) = - 



e lS < sin 5i 
Mp 



(3.17) 



The 1 So phase shifts calculated in this way is shown in the upper left panel of Fig. 01 where they are compared to 
those found by fitting Cs, as described in the previous section. The agreement between the two methods is very 
good, which confirms the ability of the subtraction to reproduce the results obtained by fitting C. The phase shifts 
produced by the CD-Bonn potential [28| are also shown in Fig. [4] Since this potential reproduces the NN data in 
this region with x 2 /d.o.f ss 1 its Si g can be regarded as a parameterization of experiments. But the CD-Bonn values 
for Si s are not well reproduced by the LO xPT potential. In this channel the LO calculation obviously has sizable 
higher-order corrections. In fact the description is known to improve when additional terms in the chiral expansion 
of V are included HHHI. 

Meanwhile the other three panels of Fig.Ushow results for the triplet channel, where we adopt the Stapp convention 
[30| for the phase shifts and the mixing parameter 



sCs^ 



arctan 



arctan 



— arctan 
2 



-z(£q 2 + tgp) 
3[(4)/co S (2e)] 



K[(i* )/cos(2e)] 
i[(t| 2 )/cos(2e)] 



3?[(4 2 )/cos(2e)] 



(3.18) 



The agreement between the two methods employed to perform the LO xPT calculation of these quantities is satisfac- 
tory, and for S3g 1 and S3£ >1 the agreement of the LO calculation with the CD-Bonn phases is strikingly good. 

The main numerical challenge of our subtractive renormalization scheme originates in solving Eq. (|3.8[) and 
Eq. (|3.16|) . Matrix elements with large p' are small and in the subtraction scheme they are calculated from a 
combination of larger matrix elements. We now compare the efficiency of the subtraction and fitting methods. We 
first note that due to the multi-step nature of subtraction, it does take about 3 times longer to run the code than 
it does for the fitting method with the same number of mesh points. On the other hand, as Table IIIII shows, the 
subtraction method converges much faster than the fitting method with respect to the number of mesh points used 
in the calculation. With 20 mesh points we already have 4 significant-figure accurancy for the subtraction method, 
while to reach this accurancy it takes more than 200 mesh points for the fitting method. The subtraction method 
converges so much faster because in the fitting method Ct affects only the part of the off-shell amplitude near p = A, 
and the oscillation phase it sets there must be communicated down to p = to get a particular scattering length a t . 
In the subtraction method the correct value of the scattering length is enforced in the equation itself, and does not 



6 



have to be achieved by fine tuning the constant Ct so as to obtain the correct behavior of the half-off-shcll amplitude 
for p of order A. 

For the same number of mesh points, converged phase shifts calculated by the two different methods differ by 0.1- 
1.3% depending on the energy, as is shown in Fig. [5] The mixing parameter for the triplet channel is more sensitive 
to numerical errors, and in that quantity there is a 1—2% relative difference between the two methods. However, we 
still claim that both methods are equivalent. 

Frederico et al. also follow the steps around Eqs. (|3.1ip and (|3.12p . although there the difference of t-matrices 

is constructed between a given cm. energy E of the NN system and a large negative energy — jj 1231 . (See also 
Refs. [3l|,[32j, where similar "first-resolvent" methods are employed.) Taking /i = in Eq. (8) of Ref. [23j| yields our 
Eq. (|3.16p . The difference between that work and what we have done here is that in Ref. [23j the Born approximation 

t = V with V given by Eq. (|2.ip is used to determine t(—^j). The final results are then independent of the "subtraction 
point" jU, provided that /x is large enough. This procedure leads to numerical results for phase shifts which are similar 
to ours. But it relies on the validity of the Born approximation at large negative energies. Since the tensor potential 
that operates in the triplet channel ~ -j at short distances the Born approximation is not actually valid at any 
energy [33l.l34j]. 

We can see this failure of the Born approximation in the triplet channel by examining the off-shell t-matrix at 

2 

a number of negative energies E = — ^ and comparing the results to V itself. Fig. [5] shows that the behavior of 

2 

tt(p,p 1 ; — jj) at fixed p' is completely different from the behavior of V. (In contrast, in the singlet channel the Born 
approximation appears to work quite well if /j, is large enough.) This defect in Frederico et a/.'s argument is manifested 
in the fact that fi dependence in physical quantities disappears only slowly as /i — > oo. The mixing parameter e±, which 
is particularly sensitive to the tensor potential, is one significant example of this. Our subtraction method gives quite 
good convergence of t\ with respect to cutoff, as shown in Fig. [7l and the convergence of other physical quantities 
with respect to A as A — > oo is quite rapid, see Tables H1HT1 (c.f. Table 1 of Ref. [23|]). This presumably results from 
the fact that our subtractive renormalization of the LS equation is based only on Hermiticity of the potential and 
momentum-independence of the contact interaction C . Under these assumptions all the results of this section and 
Sec. IIVI below are determined by the form of one-pion exchange and the value of the NN scattering lengths that are 
supplied in Eqs. and (|3~5|) . 



IV. SUBTRACTIVE RENORMALIZATION IN BOUND-STATE CALCULATIONS 

The neutron and proton form a bound state in the 3 Si- 3 Di partial-wave state: the deuteron. In this section we 
review the methods by which the deuteron binding energy is deduced from the NN t-matrix, and employ the t-matrix 
obtained via subtractive renormalization in the previous section to obtain NN binding energies and wave functions. 

Let us consider the Schrodingcr equation 

Hi> = Eij) = -Bij) (4.1) 

2 

with H = Hq + V, Ho = the free Hamiltonian, V the previously rcnormalizcd potential, and B the unknown 
binding energy. By rearranging and partial-wave expanding Eq. (|4.1[) and defining: 

(p|r,) = (-b - 4)(pI^> = (p\(- B - HoM) = (plg'H-BM), (4.2) 



M ' 



we obtain 



< 



p|r ; >=- f dp'p' 2 v u ,(p,p') 1 — JT <p'\T v > . (4.3) 



The kernel vw(p,p') {—B — ^j)" 1 of Eq. (|4.3p is the same as the kernel of the equation for the t-matrix (|2.3p . We 
see from Eq. (|4.3[) that the kernel will have an eigenvalue X(E) = 1 when the energy is equal to that at which a 
bound-state occurs, E = —B. Thus by varying B and solving the eigenvalue problem we can obtain the binding 
energy If the constant Ct is fit to a scattering length of at = 5.43 fm we obtain B = 2.13 MeV for A = 50000 McV. 
In Table ITVl the values of B = Bfu for different cutoff parameters are listed. We note that for small cutoff parameters 
A, the fitting procedure leads to a slightly larger value of Bf it . In Fig. [8] the trend of eigenvalues A e is shown as a 
function of Bf u for A = 50000 MeV, and we see that the deuteron is actually the 8th excited state of the system for 
this cutoff. As A — > oo the number of bound states in the NN system increases monotonically, but renormalization 
ensures that the shallowest one is always at B w 2 MeV. 
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We now want to demonstrate that we can apply the approach discussed in Sec. 3 to the bound-state problem. 
Multiplying Eq. (|3TTB|) with (E + B) we obtain 



(E + B) t{E) = i(0) (E + B) + (g(E) - g(0)) (E + B) t{E) 



(4.4) 



When considering the bound state, i.e. 
since, the t-matrix has a pole at E = 



E 



> — B. the first term in the brackets on the right-hand side vanishes. But, 
we can define the residue at the pole as f(—B) = lim [E + B) t{E). 

E — > — B 

Consequently, upon taking the limit, Eq. (|4.4|) becomes 



-B, 



f(-B) = i(0) g(-B) - ,g(0) f(—B) = i(0) 



1 



1 



-H n -B 



-H 



f(-B). 



(4.5) 



If this equation is written for the partial-wave form of the t-matrix then the residue operator fu'(—B) is related to 
the function (p\Ti) of Eq. flOJ) by 



( P \r l )(r l ,\p') = (p\f w (-B)\ P '). 



(4.6) 



We already calculated i(0) in Section 3, thus the kernel of Eq. (|4.5[) can be evaluated directly. The energy E = —B 
at which it has an eigenvalue of 1 reveals the bound-state energy, and the eigenvector is then the function (p\Ti) 
(up to an overall constant). Our numerical calculation gives a deuteron binding energy B = B su b =2.11 MeV. This 
result is independent of the cutoff to better than 0.5% for A > 500 MeV (see Table ITVf . The relative difference in 
the binding energy of the shallowest state obtained by solving Eq. (|4.3p and Eq. (|4.5p is 0.94% for A = 50000 MeV. 
From this we conclude that the subtractivc renormalization works for the bound state as well (or better) than it does 
for scattering observables. Fig. [5] shows the trend of eigenvalues A e of the kernel of Eq. (|4.5p as a function of B su i,. 



Our calculation gives a deuteron binding energy different from the experimental value -B e xp=2.2246 MeV. The reason 
for this difference is that we work only with the LO xPT potential, and so the analytic continuation of the scattering 

amplitude from the scattering to the bound-state regime is only accurate up to corrections of order " 2 . 

In spite of this discrepancy with experiment we can still check that our eigencnergy is the correct one for LO xPT 
with the experimental value of at reproduced. To do this we use the S-wave effective-range expansion for the on-shell 
t-matrix tg , 



m4 = 



777 + hnpl-ipo 



(4.7) 



with r t being the triplet effective range. This form of £q is valid for on-shell momenta po that obey \po\ < The 
pole in the t-matrix representing the shallowest bound state of the NN system is inside this radius of convergence. 
Thus we can solve 



1 1 2 

o r *Po + l Po = 

a t i 



(4.8) 



for po = i^f in order to determine its location, and B is then equal to j 2 /M. Using the experimental value of at, and 
the value of r t extracted from our phase shifts, we obtain B = 2.07 MeV. This value has a relative difference of ~ 1% 
with respect to the two previously calculated values. This is consistent with the omission of terms of 0(pq) in the 
effective-range expansion form employed in Eq. (|4.7[) . We now turn to extraction of the deuteron wave function. In a 
particular partial wave this is related to the above defined function (p\Ti) via 



iplip) 



B 



(4.9) 



Since (p\Ti) is equal to the eigenvector, up to an overall constant, it is now a simple matter to obtain tpi{p). In 
the upper part of Fig. [TO] the momentum-space S-wave and D-wave calculated from the above subtraction method 
are compared with those obtained from the CD-Bonn potential. The overall constant is fixed by our normalization 
condition 



dp P 2 [^{p) + ^(p)} = 1. 



(4.10) 



Here ipoifo) is the 3 S'i( 3 Di) wave function. The agreement of the LO X-PT momentum-space deuteron wave function 
with that obtained from the CD-Bonn potential is remarkable, especially in the S-wave. 
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Finally wc obtain the coordinate space wave functions by Fourier transformation 



r 




dpp 2 jL{pr)ipL{p), 



(4.11) 



where jh is the spherical Bessel function, with L = for the 3 Si wave and L = 2 for the 3 Di wave. In bottom 
part of Fig. [TO] the wave functions u(r) = uo(r) and w(r) = U2(r) are shown as a function of r. We see that our 
coordinate-space wave functions have a node near the origin, which is due to the fact that the LO X-PT potential 
has additional bound states at larger values of B, i.e. the B = 2.1 MeV state is not the ground state. The short- 
distance behavior of the deuteron wave function obtained from the LO xPT potential has been derived analytically, 
and discussed in considerable detail, in Ref. [35| . Our results for u(r) and w(r) at the distances where the Fourier 
transform is numerically stable appear to be in agreement with that work. 

In fact, Figs. [H] and [5] show that for this value of the cutoff another eigenvalue crosses 1 at B = 1.76 GeV. This is the 
binding energy of the second-shallowest bound state for A = 50000 MeV; B for this state is slightly smaller at lower 
cutoffs. The binding momentum of this state is about 1.2 GeV, and so this state is outside the radius of convergence 
of xPT: the theory's prediction of its existence is not reliable, and should be disregarded. For completeness, we show 
the wave function corresponding to this binding energy in Fig. 1111 The momentum distribution confirms that the 
prediction of this state by the low-momentum effective theory is not reliable. 



We have developed a subtractive renormalization of the S-wave Lippmann-Schwinger equation for the leading-order 
NN potential of xPT. Our results are equivalent to earlier work where the contact interaction in NN S-waves was 
fitted to obtain a particular scattering length. The advantage of our subtraction formulation is that all reference to 
the unphysical, regularization- and rcnormalization-dcpendcnt contact interaction C disappears from the scattering 
equation. Instead of searching for the unknown constants Cs and Ct, we directly use the NN scattering lengths a s 
and at as input and thereby obtain low-energy phase shifts as well as the spectrum and wave functions of all NN bound 
states which are within the domain of validity of xPT. Finding Cs and Ct for large cutoffs A can be numerically 
challenging. Subtractive renormalization avoids this problem, and so allows us to compute NN phase shifts from the 
LO xPT potential up to A = 50 GeV — a factor of 10 larger than previously employed in momentum-space calculations. 

The subtractive technique hence yields phase shifts — and bound-state information — based only on four assumptions: 

1. The form of the long-range potential (in this case assumed to be one-pion exchange); 

2. The validity of the Lippmann-Schwinger equation; 

3. Hermiticity of the underlying NN potential; 

4. Energy independence of the short-distance piece of the NN potential — or at least approximate energy dependence 
over the range of energies considered. 

Our results therefore rest on fewer assumptions than were employed to get NN phase shifts in the subtractive approach 
of Ref. [23j . In particular, we do not assume that the Born approximation is valid for the unregulated NN potential 



We have chosen to develop a subtraction based on using the NN scattering lengths as input. In principle experimental 
data at any energy for which the effective theory is valid could be used to determine T(pq,pq; E) (j>q — ME), however, 
E = seems a natural choice since it introduces the NN low-energy parameters directly into the LS equations for the 
s-waves. This freedom to choose the subtraction point was exploited to develop Callan-Symanzik-type renormalization- 
group equations for the NN system in Refs. [H, [24| . 

For higher partial waves such as the NN P-waves ( 3 Po, 3 Pi, 3 P2- 3 F2), we expect that similar logic will allow us to 
obtain subtractive equations there too. Furthermore, if we include the two-pion exchange potential [6|, 0] , which is the 
dominant correction to one-pion exchange in the xPT expansion for the NN potential, then we might hope that our 
results for phase shifts will be valid to higher energies. The subtraction technique will still allow us to rcnormalizc 
the scattering equation for this long-range potential, as long as the assumption that the constant C is approximately 
energy independent is still valid over the energy range considered (see also Refs. [H, HU). The ability of subtractive 
renormalization to facilitate higher-order calculations in this way has already been demonstrated in the three-body 
problem in an effective field theory with contact interactions alone [37], [H| . 



V. 



SUMMARY AND CONCLUSIONS 
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A [MeV] 


r s [fm] 


r t [fm] 


500 


1.556 


1.582 


1000 


1.441 


1.614 


5000 


1.365 


1.596 


10000 


1.356 


1.595 


50000 


1.349 


1.594 



TABLE I: Effective range r s (singlet) and rt (triplet) obtained using the method of Section 2 for various cutoffs A. 



A [MeV] 


C s [MeV" 2 ] 


C T [MeV" 2 ] 


500 


-6.069 x 10" 


e 


-2.307 x 10" 6 


1000 


-4.926 x 10" 


6 


-1.838 x 10" 5 


5000 


-3.902 x 10" 


6 


1.310 x 10~ 6 


10000 


-3.753 x 10" 


6 


2.040 x io~ 8 


50000 


-3.627 x 10" 


6 


-4.201 x 10 -8 



TABLE II: The constants Cs and Ct as determined by fitting the scattering lengths a s and at for different cutoff parameters 
A. 



Mesh points 


5( 3 Si) Fitting 


5( 3 Si) Subtraction 


20 


139.2453 


104.6811 


30 


113.8653 


104.6649 


40 


107.9943 


104.6622 


50 


105.9794 


104.6618 


100 


104.2317 


104.6616 


150 


104.0749 


104.6616 



TABLE III: The convergence of the phase shift <5( 3 Si) at laboratory kinetic energy Ti a b=10 MeV as a function of the number 
of Gauss-Legendre quadratures for the fitting and subtraction methods. Note that the fitting method converges to 104.1 at 
about 200 mesh points, so the final results of these two methods differ by w 0.5% as per Fig. 5. 



A [MeV] 


B fit [MeV] 


B suh [MeV] 


500 


2.10 


2.11 


1000 


2.12 


2.11 


5000 


2.11 


2.11 


10000 


2.11 


2.11 


50000 


2.13 


2.11 



TABLE IV: The binding energy B for different cutoff parameters A. The binding energy Bfu is obtained from the fitting 
method described in Section 2, and the binding energy B 3U b is obtained from the subtractive method given in Section 3. 
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T.JMeV] TJMeV] 

FIG. 1: (Color online) The NN phase shift d^So) and <5( 3 Si) as function of the laboratory kinetic energy Ti ab < 80 MeV for 
cutoff parameters A ranging from 0.5-50 GeV. The upper left panel shows the phase shift with Cs = 0, while the upper right 
panel shows the value with Cs adjusted at each value of A to reproduce the experimental value of a a . The bottom left panel 
shows the phase shift with Ct = 0, while the bottom right panel shows the value with Ct adjusted at each value of A to 
reproduce the experimental value of at. Note that the phase shift in both lower panels is given modulo n only, since multiple 
deeply bound states are present. 




A [MeV] 



FIG. 2: The renormalization constant Ct for the triplet channel as a function of the cutoff A, where A ranges from 0.5-4 GeV. 
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FIG. 3: (Color online) The comparison of the two renormalization methods for zero energy off-shell matrix element tu* (p,p'~, 0) 
in the triplet channel as a function of the momentum p for fixed p' — 100 MeV. 



Fitting renormalization 

— Subtractive renormalization 
■■■■ CD-Bonn 





FIG. 4: (Color online) The comparison of two renormalization methods for the lowest NN singlet and triplet phase shifts and 
mixing parameter as a function of the laboratory kinetic energy TLb < 80 MeV. Here A = 50 GeV is used. The phase shifts 
obtained from the CD-Bonn potential are also shown (dashed lines). 
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0.005 




FIG. 5: The relative difference of NN 5( 1 So) and <5( 3 Si) phase shift between two renormalization methods as a function of the 
laboratory kinetic energy TLb < 80 MeV. Here A = 50 GeV is used, and 100 Gauss-Legendre quadrature points are chosen for 
the solution of the LS equation. 



Singlet 



Triplet 



-5.0x10 




-4.0x10" 



v(p,100) 
|l=800 MeV 
|X=2000 MeV 
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|X= 12000 MeV 
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p [MeV] 




2000 4000 6000 8000 10000 
p [MeV] 



FIG. 6: The t-matrix elements t s (p, 100; —37) (left) and tg (p, 100; —^7) (right) as functions of the half-shell momentum p 
compared to v a (p, 100) and vg (p, 100). The parameter fi ranges from 0.8-12 GeV. 
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FIG. 7: The mixing parameter ei obtained with the subtraction method as a function of the laboratory kinetic energy Ti a b < 
20 MeV and for cutoff parameters A ranging from 0.5-50 GeV. 
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FIG. 8: (Color online) The eigenvalues A e obtained from Eq. (|4.3(l as a function of the binding energy B. Note that the vertical 
axis is on a log scale. 
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FIG. 9: (Color online) The eigenvalues A e obtained from Eq. (|4.5[1 as a function of the binding energy B. 
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FIG. 10: (Color online) Upper panel shows the momentum-space wave functions of the shallowest AW bound state in the 
triplet channel as a function of p, where V'o(p) is the 3 Si wave function and ^(p) denotes the 3 Di wave (solid lines). These 
wave functions are obtained from the subtracted integral equation with A = 50 GeV. Bottom panel shows the coordinate-space 
wave functions as functions of r, where u(r) denotes the 3 Si wave and w(r) denotes the 3 Di wave (solid lines). The dotted 
lines indicate the corresponding wave functions obtained from the CD-Bonn potential. 
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p [MeV] p [MeV] 



FIG. 11: Wave functions of the second-shallowest NN bound state in the triplet channel as functions of p, where i/>o(p) denotes 
the 3 Si wave and tfeCp) denotes the 3 Di wave. These wave functions are obtained from the subtracted integral equation with 
A = 50 GeV. 



